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. We have numerically investigated statistical properties of the so-called interoccurrence time or the 

f"^ ■ waiting time, i.e., the time interval between successive earthquakes, based on the two-dimensional 

(2-D) spring-block (Burridge-Knopoff) model, selecting the velocity-weakening property as the con- 
CNJ , stitutive friction law. The statistical properties of frequency distribution and the cumulative dis- 

tribution of the interoccurrence time are discussed by tuning the dynamical parameters, namely, a 
stiffness and frictional property of a fault. We optimize these model parameters to reproduce the 
interoccurrence time statistics in nature; the frequency and cumulative distribution can be described 
by the power law and Zipf-Mandelbrot type power law, respectively. In an optimal case, the h- value 
f"^ ' of the Gutenberg- Richter law and the ratio of wave propagation velocity are in agreement with those 

£f} , derived from real earthquakes. As the threshold of magnitude is increased, the interoccurrence time 

distribution tends to follow an exponential distribution. Hence it is suggested that a temporal se- 
quence of earthquakes, aside from small-magnitude events, is a Poisson process, which is observed 
| in nature. We found that the interoccurrence time statistics derived from the 2-D BK (original) 

model can efficiently reproduce that of real earthquakes, so that the model can be recognized as a 
realistic one in view of interoccurrence time statistics. 
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I. INTRODUCTION 



e 
e 

Earthquakes are complex phenomena, involving relative motion of faults. Many fundamental problems remain 
elusive, such as the source mechanism, a physical background, and so forth. To elucidate the source mechanism of 
O . earthquakes, it is essential to ascertain the physical background of the friction force acting on the surface of a fault. 
However, the unified friction law of a fault has not yet been established. Nowadays, some constitutive friction laws 
have been proposed based on laboratory experiments [l|, e.g., the velocity- weakening property 0, H|, the rate and 
O^l ' state dependent friction law [3 , and the slip dependent constitutive law [|| . On the other hand, there are well-known 
and relevant empirical laws Q for earthquakes. The most familiar one is the Gutenberg-Richter (GR) law which 
describes a relation between the seismic magnitude M and its frequency n as n oc 10 -bM , where b stands for the 
£^ | &-value and is similar to unity. 

. The time interval between successive earthquakes is often called the interoccurrence time or waiting time. Very 
recently, statistical properties of the interoccurrence time have been studied [1, 0, [lj| [HI, O EH, EH • In Southern 
California, the probability density of the return time for a given xy-region of size L follows a power law Q, which is 
governed by three factors 0: b- value of the GR law, the Omori law for aftershocks [TB|, and the fractal dimension 
\ of faults. Moreover, interoccurrence time distributions for large earthquakes exhibit both a gamma distribution and 
an exponential distribution [icl |. Corral showed after examining a wide variety of global earthquake catalogs that 
, | the probability density functions obey a generalized gamma distributi on llll . Abe and Suzuki found that cumulative 
distributions are governed by the Zipf-Mandelbrot type power law [lfj, which is equivalent to the (/-exponential 
distribution with q > 1 [l2j] based on the nonextensive statistical mechanics proposed by Tsallis [lj]]- Recently, 
Saichcv and Sornette demonstrated that the distribution of interoccurrence time is derived from the known laws of 
the GR law and the Omori law [I103. 

In the 1980s, Bak et al. proposed the concept of self-organized criticality (SOC) [IH, [l9[, according to which a 
nonequilibrium open system gradually evolves into a critical state, where a distribution of physical quantity follows the 
power law. A crust is a nonequilibrium open system because energy is supplied by plate motion and dissipated by an 
earthquake. Introducing a sand-pile model involving SOC and using the cellular automaton (CA) simulation method, 
they demonstrated an empirical power law which is related to the GR law [2fi| . They suggested that earthquakes can 
be categorized into SOC phenomena. Inspired by their study, a number of SOC-based earthquake models have been 
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proposed to reproduce statistical properties of earthquakes [21|. For example, a spring-block [Burridgc-Knopoff (BK)] 
model [13] has been widely utilized. Based on the one-dimensional (1-D) BK model with nonlinear friction force, 
hereafter denoted by BK (original) model, Carlson and Langer showed the power law distribution like GR law, in a 
restricted mag nitude region [23l . |24| . It has been shown that the power law distribution only concerns small cree pin g 
events [25], |26| Carlson and Langer also discussed the interoccurrence time statistics for large earthquake events \2l\ . 
Based on the 2-D BK (original) model extended by Carlson, the magnitude distribution [28[ and statistical properties 
of the shear stress drops [29[ were reported. This model calculation has also been performed by the CA simulation 
method, from here on referred to as the BK (CA) model, and the statistical prop erties of earthquakes, such as the 
magnitude distribution [10, HH, HI, HH and the interoccurrence time statistics [34 |35| were discussed. Note that the 
1-D BK (CA) model is not a good model for the Gutenberg-Richter law, see, for instan ce, 13611 ■ The studies of the 
interoccurrence time statistics based on its modified versions were reported \37\ , l38l . [3^ . |40L l4l| . The fundamental 
difference between the BK (original) model and the BK (CA) model is how to dissipate the friction energy; the 
friction force is governed by the constitutive law for the BK (original) model, and by the equipartition law for the BK 
(CA) model. The BK (original) model still has the great advantage of involving a constitutive friction law that well 
accounts for the laboratory experiment as well as the inertia. 

The BK model in itself traces back to the late 1980s. Wc realize that the model simplifies the complex fault 
dynamics, so that some features of the fault, such as no radiation damping and long-range interactions, arc inherently 
neglected. However, since the model can efficiently extract the statistical properties of earthquakes, such as the GR 
law and the Omori law [421] . it has been attracted much attention. Recent reports based on the BK (original) model 
have focused on the long range stress transfer [IH , a fractal structure of faults 0, [45| , the statistical properties of 
the magnitude distribution and the interoccurrence time based on the rate and state dependent constitutive law [4(| , 
the correlation of seismicity [4?], [H, [4^, and the epicenter distance statistics [50]. However, the comprehensive 
application of the 2-D BK (original) model to the statistical properties of the interoccurrence time still remains to be 
done. According to the recent papers [ID, [3] , the interoccurrence time statistics derived from the modified versions 
of the BK model [10, HO, EH, [42j which show the GR law and the Omori law reproduce these statistics in nature. 
On the other hand, in the case of the 2-D BK (original) model, whether the natural interoccurrence time statistics 
can be extracted or not is still contoroversial, because this model cannot produce aftershocks. In this study, we have 
attempted to work out how the interoccurrence time statistics are influenced by the variation of the major physical 
quantities, such as stiffness and frictional parameters of faults, and a threshold value of magnitude and by seismicity 
without aftershocks. Here, we report numerical investigations on the interoccurrence time statistics based on the 2-D 
BK (original) model by testing various dynamical parameters. These physical parameters are restricted or optimized 
so as to reproduce the statistics of actual earthquakes. Wc also compare the statistical properties derived from this 
model to those from the BK (CA) model. We show that the 2-D BK (original) model can successfully reproduce the 
interoccurrence time statics in nature. 



II. MODEL 



In Fig. 1 (a), we schematically illustrate the 2-D spring-block model which we simulated in this study. The model 
is composed of blocks, two plates, two different kinds of coil springs, and k%, and a leaf spring, k p , corresponding 
respectively, to segments of a fault, geological plates, a compression stress, and a shear stress. The fundamental idea 
of the model is that earthquakes are caused by the stick-slip motion of the fault. It is noted that this model is virtually 
the same as Otsuka's model [HI but is different from Carlson's (fcjf = kf!) [H|. We assume that the slip direction of 
blocks is restricted only to the y-direction. The equation of motion at site can be written as 



Vi 



dt 2 



Ki.Vi+1,0 + Ui-1,3 - 2 2/ij) + k c(yi,3-l + Vij+1 ~ tyhj) - k pVi 



F 



dyi,j 
dt 



(1) 



where m and j/jj are, respectively, mass and displacement of the block and F is a dynamical force between the block 
and the bottom plate. 

Now we rewrite Eq. (1) into a dimensionless form. The dimensionless dynamical friction force (f> is obtained as 
4>{y/ v i) = F(y)/FQi where Fq and v\ are the maximum friction force and the characteristic velocity, respectively. A 
dimensionless time t' and displacement Ui.j arc defined by 



t' = ui p t = J kp/m t, Uij = = 



Vi,3 
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Then the dimensionless form is 
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FIG. 1: (a) 2-D spring-block model. In our calculation, the system size is (100 x 25) on the (x, y) plane, fcj, k]i, and k p are 



spring constants. The friction force acts on the surface between the block and the bottom plate, 
nonlinear dynamical friction function, a = 0.01 throughout the simulation. 



(b) a-dependence of the 



and. 
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where v is the slipping velocity, being on the order of 1 m/s in a real crust. In this work, we choose the friction force 
as the velocity-weakening constitutive law which states that as the slipping velocity increases, the dynamical friction 
force decreases. This constitutive law was observed in the accurate rock fracture experiment [52j . Then, we use the 
dimcnsionlcss dynamical friction force <p given by 

f (-oo,l], L> = 0, 

4,(11) = } (J>Q (3) 

I {l + 2a[f7/(l-<7)]}' 

The formulation given in Eq. (3) was introduced by Carlson et al. [25[ [see Fig. 1 (b)]. This can be regarded as an 
ideal friction function involving the stick-slip motion, so that it has often been adopted to 1-D and 2-D BK (original) 
models 0, [H, S3, 0, H, [H, [13, 0, [H, [H3l ■ To exclude a back slip, (U tJ < 0), we treat 4> as a tunable parameter 
when U = 0. 

The model is governed by five parameters, l x ,l y ,a,a, and v. l x and l y are dimcnsionlcss stiffness parameters in 
the x- and y- direction, respectively, a means a decrement of the dynamical friction force with increasing slipping 
velocity, a is a stress gap between the normalized maximum friction force (= 1) and dynamical friction force 0(0). v 
is a dimensionless loading velocity. l x and l y are related to Lame's constants, A and fj,, as [5l| 

«-(£)'■«-££(£)'■ 

where Ax, Ay, and Az are infinitesimal lengths in the x-,y-, and z- directions, respectively. According to the experi- 
ments P, [52|, a is positive and on the order of 10° (=1). We can treat v as being equal to zero, when a slip event 
occurs because v, in nature ~ 10~ 9 , is far smaller than a slipping velocity of ~ 1. This approximation ensures that 
no other event takes place during an ongoing event. Obviously, v is related to loading time and interoccurrence time. 
We fixed v = 0.01 because we will rescale it later. 

In this paper, Ax = Ay = Az is used with assumptions of A = /i. We solved Eqs. (2) and (3) by the fourth-order 
Runge-Kutta method under a free boundary condition. A small irregularity of block displacements is considered at 
an initial stage of the calculation. We use a 10 5 order of earthquakelike events after some period when the initial 
randomness does not influence the statistical properties. Then the interoccurrence time statistics arc systematically 
studied by changing l x and l y , and by a increasing the threshold of magnitude M c . 
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FIG. 2: The frequency distributions of the interoccurrence time for different a (A : a — 1.5, □ : a = 2.5, Q : a = 3.5, and 
+ : a — 4.5). f is the characteristic interoccurrence time, with fixing l x — 1 and l y — \/3. b' is calculated from the slope of the 
dashed line. All plots except for in the case of a = 4.5 are shifted vertically for clarity. 



III. RESULTS AND DISCUSSIONS 



In this model, a slip of blocks is considered as an earthquake. An event starts when a block begins to slip in the 
direction toward +y and ends when all blocks stop slipping. We define the interoccurrence time as the time interval 
between successive events. Accordingly, nth interoccurrence time can be described as r n = t' n — tf n _i, where t' n and 
if n _x are nth and n — 1th earthquake occurrence time, respectively. 



A. Friction parameter a dependence 



In the first performance of our simulations, we study the a- dependence of the interoccurrence time statistics. For 
that purpose, the stiffness parameters l x and l y are fixed (l x = 1 and l y = v3), and the friction parameter a is changed 
from 1.5 to 4.5. As a result, the mean value of the interoccurrence time, T, gradually increases as a is increased; 
namely, T ~ 1.73 (a = 1.5), 2.65 (a = 2.5), 3.40 (a = 3.5), and 3.53 (a = 4.5). The frequency distributions of 
the interoccurrence time for different a are shown in Fig. 2 as a function of a scaled interoccurrence time, r' = r/f, 
where f represents characteristic interoccurrence time and is set at 2.0 arbitrarily. We confirmed that in the case of 
a > 4.5, the statistical properties do not change qualitatively. The frequency distributions of any a exhibit a power 
law in the short-time region of 0.5 J$ r' ;$ 7. The power law exponent b' gradually decreases as a is increased, for 
example, b' ~ 3.04 (a = 1.5), b' ~~2.40~(a = 2.5), b' ~ 1.96 (a = 3.5), and b' ~ 1.78 (a = 4.5). As shown in 
Fig. 2, the distributions are categorized into three types; first, in the case of a = 1.5 and 2.5, the frequency of a long 
interoccurrence time region of r' > 10 is enhanced more than expected by the power law decay, hereafter denoted 
by type A. Second, for a = 4.5, the frequency becomes less than predicted by the power law, which is referred to as 
type B. Finally, for the intermediate case of a = 3.5, the distribution most closely follows the power law, which is 
described as type C. It is suggested that the system undergoes a critical state when a ~ 3.5 because the distribution 
shows the power law. We note that a critical power law has never been achieved qualitatively using the BK (original) 
model, whereas it has generally been obtained with the 2-D BK (CA) model [4fj| . 

Now, we discuss statistical properties of cumulative distributions of the interoccurrence time upon changing a. 
According to (l2j |. the cumulative distributions of the interoccurrence time, denoted by P{> t'), with earthquake 
data follow the Zipf-Mandelbrot type power law, 

P{> r') = * = M-t'/to) = [(1 + (1 - 9 )(-t'/to))^] + , (5) 
(1 + er)~ 

where ([a]+ = max[0, a]) and £,7,9, and to are positive constants. Especially, q, tq and e q (x) are called an entropy 
index, a time-scale parameter, and the ^-exponential function, respectively. e q (x) converges to an exponential function, 
e x , as q — ► 1. Figures. 3 (a) and 3 (b) are examples of P(> t') for a = 2.5 and 3.5. In order to estimate an optimal a, 
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FIG. 3: The cumulative distributions of the interoccurrence time for a — 2.5 and 3.5 with fixing l x = 1 and l y = a/3- In this 
figure, the dotted, solid, and broken lines correspond, respectively, to the simulation data, the Zipf-Mandelbrot type power law 
defined in Eq. ((5]), and the exponential distribution. 



where the model yields or reproduces the interoccurrence time statistics in nature, we fit our numerical data to the Zipf- 
Mandelbrot power law given by Eq. (5). As a result, these fitting parameters are estimated to be q = 1.10, tq = 2.26, 
and p z = 0.986 for a = 2.5 (a) and q = 1.08, r = 2.05, and p z = 0.989 for a = 3.5 (b) by the least-squares method. 
p is a correlation coefficient between our simulation data and the fitting curve. For both a = 2.5 and 3.5, q turns 
out to be closed to unity, so we introduce a different fitting model of the exponential distribution function described 
by P(> t') = Aer T l T . The fitting parameters evaluated are A = 0.16, f = 5.08, and p e = 0.988 for a = 2.5, and 
A = 0.18, f = 3.97, and p e = 0.980 for a = 3.5. The data demonstrate that in the case of l x = 1 and l y = s/3, 
P{> t') can be best described by the Zipf-Mandelbrot type power law, simultaneously exhibiting quantitatively the 
nature of the real seismicity. Additionally, P(> r') gives a good description of the Zipf-Mandelbrot type power law 
in the case of a small fe'-value, judging from the small difference between our calculation and an ideal distribution, 
around 0.5 J$ r' ;$ 3. We note that the time-scale parameters, To and f , are in particular sensitive to the variation of 
a, compared to the other fitting parameters. 



B. Stiffness parameter l x and l y dependence 

In the second run of our simulations, the stiffness parameter-dependence of the interoccurrence time statistics is 
investigated, whereas a is fixed at 3.5, and l x and l y are systematically changed. We study the two different cases. 
One is an isotropic case, l x = l y , and the other is an anisotropic case, such as l x = 1 and l y = v3- As a result, 
T depends on l x and l y , for instance, f ~ 3.94 (l x = l y = 1), 2.24 (l x = l y = 0.973 (l x = s/3 and l y = 3), 

and 2.08 (l x = y/2 and l y = yS). In Fig. 4, we display the frequency distributions for different stiffness parameters. 
The 6'-values also depend on the stiffness parameters, such as b' ~ 2.67 (l x = l y = V3), b' ~ 3.09 (l x = V3 and 
l y = 3), b' ~ 2.65 (l x = y/2 and l y = Vb), and b' ~ 1.94 (l x = 1 and l y = \/T0). Thus we point out that among the 
systematically varied stiffness parameter settings, in almost all cases, the frequency distributions are categorized into 
type A, exhibiting a board peak in the long time region of 10 $ r' $ 30. This implies that quasiperiodic events occur, 
which has been reported in the 1-D BK (original) model 2^, 4?], [48j and in the 2-D BK (original) model [2|| , except 



for the case of l x = l v = 1. 

As shown in Fig. 5, we discuss the cumulative distributions for different stiffness parameters, such as l x = l y = s/3 
(a) and l x = y/3 and l y — 3 (b). The Zipf-Mandelbrot type power law and the exponential distribution are used 
as distribution functions to fit our simulation results. Then, the Zipf-Mandelbrot power model yields q = 1.08 and 
t = 2.67 for (a), and q = 1.13 and To = 1.23 for (b), giving a correlation coefficient p z = 0.954 and 0.921, respectively. 
As for the exponential distribution, A = 0.08 and To = 5.52 for (a) and A = 0.017 and To = 4.46 for (b) are deduced 
together with p e = 0.974 and 0.962, respectively, for (a) and (b). 
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FIG. 4: The frequency distributions of the interoccurrence time for different l x and l y (A : l x = 1 and l v — vTO, □ : l x = \/2 
and Z y = \/5, O ■ lx = l y = V%, and + : Z^ = v^3 and l y = 3), with a = 3.5. All plots except for in the case of l x = 1 and 
ly = VW are shifted vertically for clarity. 




FIG. 5: The cumulative distributions of the interoccurrence time for l x — l y — \/3 in (a) and l x = v3 and l y — 3 in (b), while 
a is fixed to be 3.5. In this figure, the dotted, solid, and broken lines correspond, respectively, to the simulation data, the 
Zipf-Mandelbrot type power law defined in Eq. ([5]), and the exponential distribution. 



C. Threshold of magnitude M c dependence 



We study the magnitude-dependence of the interoccurrence time in order to discuss the interoccurrence time 
statistics for large magnitude earthquakes. In this model, we defined a seismic moment Mq and a seismic magnitude 
M as Mo = Y17j fiUij and M = (logMo)/1.5, where SUij is a total displacement at site during an event and n 
is the number of slipping blocks (53j . One may find it unrealistic at first sight that the seismic magnitude M can be 
negative. However, it is natural in the case where n J$ 10, because the dimensionlcss displacement per block is less 
than unity by definition. Here, we examine the interoccurrence time statistics by altering M c , using the simulation 
data for l x = l,l y = and a = 3.5. It should be noted that the upper limit of M c is optimized so as to ensure 
sufficient data for us to evaluate the statistical properties. 

For the parameter setting given above, the seismic magnitude M ranges from —0.86 to 1.5. In order to keep the 
statistical quantities accurate, we determined the upper limit of M c to be —0.60. Then, we selected —0.80, —0.75, 
—0.70, —0.65, and —0.60 as suitable M c values to test. The frequency distributions M c = —0.70 and —0.60 are 
presented in Fig. 6 (a). The frequency of short interoccurrence times satisfying 0.5 J$ t' ;$ 2 gradually decreases 
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FIG. 6: The interoccurrence time statistics increasing the threshold of magnitude M c , for l x = l,l y = v3, and a = 3.5 
(O : M c = —0.70, and A : M c — —0.60). The solid line corresponds to the exponential distribution. The frequency and 
cumulative distributions are in (a) and in (b), respectively. 



TABLE I: Summay of the interoccurrence time statistics based on the 2-D BK (original) model and on the real seismicity. 
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"As we mentioned before the frequency distributions are categorized as three types, namely, type A, type B, and type C. 
''Two test distribution functions are introduced, namely the Zipf-Mandelbrot power law and the exponential distribution denoted here 
the power law and exponential law, respectively. 



when we increase M c . This indicates that small magnitude events occur successively. It is found that the frequency 
distributions follow the exponential distribution for large M c . This implies that the cumulative distributions are 
also described by the exponential distribution. As shown in Fig. 6 (b), indeed the exponential law does govern the 
cumulative distributions for large M c . Note that based on the 2-D BK (CA) model, the distribution for large M c is 
governed by the power law [38| , which differs from our results for the 2-D BK (original) model. 

We can demonstrate that the interoccurrence time statistics depends on the threshold of magnitude M c ; the 
frequency and the cumulative distributions are definitely changed qualitatively and quantitatively as M c is increased. 
For large M c , these distributions follow the exponential distributions. Moreover, it is found that a temporal sequence 
of events, except for small magnitude events, is a Poisson process, which is also reported using real earthquake 
data [n|. We point out that the transition magnitude point can be estimated to be —0.65. 

D. Optimization of the model parameters and comparison with the real seismicity 

Now, we compare the interoccurrence time statistics obtained from our model with that observed in nature in order 
to confirm our model is realistic in view of extracting the interoccurrence time statistics of an actual earthquake. Here, 
the dynamical parameters are tuned so as to extract the statistical properties of earthquakes as clearly as possible. 
As displayed in Table I, the optimal parameters are estimated to be l x = 1, l y = -\/3, and a = 3.5. It is noted that the 
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FIG. 7: The interoccurrence time statistics obtained from the model and real earthquake data. The frequency distribution and 
the cumulative distribution are shown in (a) and (b), respectively. In (b) the fitting parameter of the cumulative distribution 
yields q — 1.13 and to = 1.72, q = 1.05 and to = 1.58, and q = 1.08 and to = 2.05 from the Southern California data, and the 
Japan data, and our simulation, respectively. 

exponent b' in nature is similar to unity because this value is related to the Omori law as b' = 2 — (1/p) [54|, where p 
characterizes the Omori law. Our findings are in agreement with real data semiquantitatively because the 6'-value is 
greater than 1. 

We display the cumulative distribution of interoccurrence time for Southern California and Japan overlapped in 
our simulation data in the case of the optimal parameters, namely, l x = l,l y = and a = 3.5 in Fig. 7 (b). It 
is remarked that in order to compare the statistical properties, we used the data from [l2| and then we reseated the 
interoccurrence time by f = 1000 s. It is found that the cumulative distributions also qualitatively satisfied those of 
real seismicity, the Zipf-Mandelbrot type power law. As can be seen from Figs. 3 (b) and Fig. 8 (b), there appears a 
small but distinguishable difference between our numerical data and the Zipf-Mandelbrot power law in the region of 
0.4 ;$ t' i5 4. When the 6'- value gradually decreases, our data of P(> r') become closed to the Zipf-Mandelbrot type 
power law. If the model reproduced the Omori law, corresponding to b' ~ 1.2, we could obtain realistic interoccurrence 
time. Therefore to obtain more realistic interoccuurence time statistics, we need to select a model which reproduces 
the Omori law as well as the GR law. 

Finally, we discuss other remaining statistical properties, such as b- value of the GR law and the ratio of the seismic 
waves. We calculate the 6-value from the slope of the power law magnitude distribution. The results are listed in 
Table. I. As for the optimal parameters, l x = 1, l y = v3, and a = 3.5, the b- value is similar to unity, which is consistent 
with the 6-value in nature. As already described, l x and l y are defined, respectively, as yjk% /k p and \Jk v c /k p . Here, 

we can rewrite them into l x = yj ^j™ aud l y = y ^"/™ - lx an d ly correspond respectively to the secondary wave 

(S-wave) and the primary wave (P-wave) because in the x- and y- direction, the oscillation direction of blocks is 
respectively perpendicular to and parallel to the wave propagation direction. The velocities of the P-wave and S-wave 
in real crust are closed to 7 and 4 km/h, respectively, so that a real observation gives l y /l x ~ 1-7. For l x = 1 and 
l y = \/3, the ratio l y /l x well coincides with that of the real earthquakes. Although for other cases, e.g., l x = \J~2 and 
l y = s/5, and l x = s/3 and l y = 3, the ratio is again similar to 1.7, the correlation coefficient is found to be worse than 
that of l x = 1 and l y = V3. Therefore we found that the 6-value of the GR law and the ratio of seismic waves satisfy 
those of earthquakes in nature in the case of optimal parameters (l x = 1, l y = \/3, and a = 3.5) which derived from 
the interoccurrence time statistics. 

E. System size dependence 

We study the size-dependence of the interoccurrence time statistics. For this purpose, the system size, N, is varied 
from 625 (25 x 25) to 22 500 (150 x 150), with fixing l x = 1, l y = V3 and a = 3.5. In this case, the model reproduces 
a realistic 6-value and the interoccurrence time statistics. 
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FIG. 8: Size dependence of the interoccurrence time statistics for the case of l x = 1, l v = y/3, and a — 3.5. We vary the system 
size N from 625 to 22 500. (A : N = 625, □ : N = 2500, O = N = 10 000, and + : N = 22 500) 

We display the interoccurrence time statistics for different system sizes, N, in Fig. 8. As clearly be seen from Fig. 
8 (a), the frequency distribution shows the power law in the region of 1 ;$ r' ;$ 20 with an exponent, b' ~ 2.07, 2.16, 
2.18, and 2.16 for N = 625, 2500, 10 000, and 22 500, respectively. This shows that the statistical property of the 
distribution, the power-law exponent b' , converges, for N 2; 2500. On the other hand, the maximum interoccurrence 
time gradually increases when TV is increased. It is shown in Fig. 8 (b) that the cumulative distributions do not 
depend on N, in the region of 0.1 ^ r' ^ 10, whereas it is found that the tail of the distributions becomes stretched 
for large N. By fitting the data to the Zipf-Mandelbrot power law given in eq. (5), the optimal fitting parameters 
are obtained; q = 1.04 and r = 2.34 for N = 625, q = 1.06 and t = 2.55 for N = 2500, q = 1.11 and r = 2.48 for 
N = 10 000, and q = 1.15 and To = 2.10 for N = 22 500. The correlation coefficient, p, between our numerical result 
and the ideal curve is found to be 0.989, 0.986, 0.989, and 0.970 for N = 625, 2500, 10 000, and 22 500, respectively. 
The result manifests that the statistical properties of P(> r') remain vertually unchanged for N ;$ 22 500. Therefore, 
we conclude that the interoccurrence time statistics hold up when 2500 iS N J$ 10 000. It should be noted that in our 
numerical simulation, the system size, N, is set at 2500, where the size-dependence is negligible. 

IV. CONCLUSION 

In light of our numerical investigations based on the 2-D BK (original) model, the interoccurrence time statistics 
depends on the dynamical parameters, stiffness, l x and l y , and frictional properties of a fault, a. Then we are able to 
restrict or optimize the dynamical parameters so as to reproduce of the interoccurrence time statistics in nature; the 
frequency distribution shows the power law and the cumulative distribution reveals the Zipf-Mandelbrot type power 
law. Additionally, we demonstrate that the interoccurrence time distributions follow the exponential distributions 
for large M c , so that a temporal sequence of earthquakes, except for small-magnitude earthquakes, is reduced to a 
Poisson process. We have found that the 2-D BK (original) model with the optimal parameters l x = l,l y = s/3, and 
a = 3.5 can be recognized as a realistic model for earthquakes in view of its fine reproduction of the major statistical 
properties of earthquakes, such as the interoccurrence time statistics, 6-value of the GR law, and the ratio of the 
seismic waves. We conclude that the 2-D BK (original) model is sufficient to extract the realistic interoccurrence 
time statistics without employing the BK (CA) model that assumes an unrealistic friction law and no inertia. The 
present results underlie the advantages of the 2-D BK (original) model for its extensive application to other mechanical 
systems involving stick-slip behavior. 
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